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This  paper  considers  compressible  turbulent  flows  at  low  turbulent  Mach  numbers.  Con¬ 
trary  to  the  general  belief  that  such  flows  are  almost  incompressible,  (i.e.  the  divergence 
of  the  velocity  field  remains  small  for  all  times),  it  is  shown  that  even  if  the  divergence  of 
the  initial  velocity  field  is  negligibly  small,  it  can  grow  rapidly  on  a  non-dimensional  time 
scale  which  is  the  inverse  of  the  fluctuating  Mach  number.  An  asymptotic  theory  which  en¬ 
ables  one  to  obtain  a  description  of  the  flow  in  terms  of  its  divergence-free  and  vorticity-free 
components  has  been  developed  to  solve  the  initial-value  problem.  As  a  result,  the  various 
types  of  low  Mach  number  turbulent  regimes  have  been  classified  with  respect  to  the  initial 
conditions.  Formulae  are  derived  that  accurately  predict  the  level  of  compressibility  after  the 
initial  transients  have  disappeared.  These  results  are  verified  by  extensive  direct  numerical 
simulations  of  isotropic  turbulence. 


'Research  was  supported  by  the  National  Aeronautics  and  Space  Administration  under  NASA  Contract 
No.  NAS1-18805  while  all  of  the  authors  were  in  residence  at  the  Institute  for  Computer  Applications  in 
Science  and  Engineering  (ICASE),  NASA  Langley  Research  Center,  Hampton,  VA  23665. 


I.  Introduction 


Turbulence  is  the  most  common  state  of  fluid  motion.  It  is  an  all  pervasive  ubiq¬ 
uitous  phenomenon  present  in,  to  cite  a  few  instances,  weather  patterns,  ocean  cur¬ 
rents,  high-temperature  plasmas,  astrophysical  jets,  and  combustion.  Despite  the 
best  theoretical  and  experimental  attempts  of  more  than  a  century,  and  more  recent 
computational  approaches,  turbulence  remains  to  be  one  of  the  great  unsolved  prob¬ 
lems  of  fundamental  physics,  and  poses  a  grand  challenge  to  any  type  of  scientific 
investigation. 

There  has  of  course  been  some  progress  in  our  understanding  of  turbulence  in 
low-speed  flows  which  are  essentially  incompressible  (see  Dwoyer,  Hussaini  and  Voigt 
1984).  There  have  also  been  a  few  attempts  at  predicting  the  effect  of  compressibility 
on  some  known  results  of  incompressible  turbulence  such  as  the  Kolmogorov  spectrum 
(Zakharov  and  Sagdeev  1970,  Kadomtsev  and  Petvishvili  1973,  Moiseev, Sagdeev,  Tur 
and  Yanovskii  1977,  and  L’vov  and  Mikhailov  1978).  The  majority  of  work  in  com¬ 
pressible  turbulence  focuses  naturally  on  the  linear  regime,  owing  to  the  inherent 
difficulty  of  treating  nonlinearities  which  now  involve  variable  density  as  well.  Moyal 
(1952)  appears  to  be  the  first  one  to  look  at  spectra  of  homogeneous  isotropic  turbu¬ 
lence  in  compressible  fluids.  He  decomposed  the  velocity  field  in  the  spectral  space 
into  a  longitudinal  component  (random  noise)  and  a  transverse  (eddy  turbulence) 
component.  This  is  in  fact  equivalent  to  a  Helmholtz  decomposition  (grad'^nt  of  a 
scalar  and  curl  of  a  vector)  in  physical  space.  The  longitudinal  component  in  physical 
space  is  variously  known  as  the  acoustic  component  or  compression  turbulence;  the 
transverse  component  is  also  termed  the  incompressible  component,  the  solenoidal 
component  or  shear  turbulence.  A  broad  conclusion  of  his  work  was  that  the  inter¬ 
action  between  these  two  components  is  exclusively  due  to  the  nonlinear  terms,  and 
such  interactions  are  the  strongest  at  high  levels  of  turbulence  and  at  high  values  of 
Reynolds  number.  The  properties  of  sound  emitted  in  the  far  field  by  eddy  turbulence 
were  studied  by  Lighthill  (1952,  1954,  1956).  His  formula  for  sound  emission  is  found 
in  laboratory  experiments  to  remain  valid  far  beyond  the  linear  regime.  Kovasz- 
nay  (1957)  proposed  a  decomposition  of  compressible  turbulence  into  three  modes  - 
the  vorticity  mode,  the  entropy  mode  and  the  acoustic  mode  -  and  showed  how  to 
determine  their  levels  and  correlations  from  mass  flow  and  stagnation  temperature 
fluctuations  measured  by  a  hot-wire  anemometer.  Chu  and  Kovasznay  (1958)  have 
outlined  a  consistent  successive  approximation  procedure  in  terms  of  an  amplitude 
parameter,  and  have  provided  explicit  formulae  for  second-order  interactions  among 
these  three  modes.  They  provide  no  explicit  solutions  since  their  main  purpose  was 
to  provide  a  consistent  framework  to  assess  the  nonlinear  interactions  in  the  experi¬ 
mental  data.  Tatsumi  and  Tokunaga  (1974)  and  Tokunaga  and  Tatsumi  (1975)  study 
the  interactions  of  weakly  nonlinear  disturbances  such  as  compression  waves,  expan¬ 
sion  waves  and  contact  discontinuities  using  a  reductive  perturbation  method  due 
to  Taniuti  and  Wei  (1968).  The  key  result  is  that  the  interaction  between  waves 


of  different  families  of  characteristics  leads  to  alterations  in  their  amplitudes,  phase 
velocities  and  propagation  directions,  whereas  the  interaction  between  waves  of  the 
same  family  of  characteristics  causes  merger  or  coalescence.  They  further  inferred  that 
the  statistical  properties  of  two-dimensional  shock  turbulence  are  similar  to  those  of 
one-dimensional  shock  turbulence  which  in  turn  are  identical  to  those  of  Burgers  tur¬ 
bulence.  More  recently,  in  a  preliminary  study  of  return  to  isotropy  in  a  compressible 
flow,  Lumley  (1989)  has  provided  some  orders  of  magnitude  estimates  in  terms  of 
Mach  number  and  Reynolds  number.  An  excellent  review  of  the  published  literature 
on  compressible  turbulence  up  to  1967  may  be  found  in  Monin  and  Yaglom  (1967). 

The  computational  approach  to  turbulence  is  based  on  the  Navier-Stokes  equa¬ 
tions.  The  first  attempt  to  solve  numerically  the  equations  of  motion  for  compressible 
homogeneous  turbulence  is  due  to  Feiereisen,  Reynolds  and  Ferziger  (1981).  They  as¬ 
sumed  the  divergence  of  the  initial  flow  field  and  its  time-derivative  to  be  both  zero. 
It  was  therefore  not  surprising  that  their  results  for  isotropic  compressible  turbulence 
did  not  show  any  departure  from  the  corresponding  incompressible  data.  Recently, 
Passot  and  Pouquet  (1986)  have  carried  out  numerical  simulations  of  two-dimensional 
homogeneous  compressible  turbulent  flows.  They  show  that  the  behavior  of  the  flow 
beyond  an  initial  turbulent  Mach  number  of  0.3  differs  sharply  from  the  lower  Mach 
number  cases  which  are  characterized  by  the  absence  of  shocks. 

In  the  present  paper,  we  develop  an  asymptotic  theory  (with  turbulent  Mach 
number  as  a  small  parameter)  to  solve  the  compressible  Navier-Stokes  equations. 
The  problem  is  set  up  as  an  initial-value  problem  to  study  the  influence  of  initial 
conditions  on  the  subsequent  turbulence  structure  and  its  dynamical  evolution.  Ex¬ 
plicit  relationships  between  the  initial  turbulent  Mach  number,  pressure  and  velocity 
fluctuation  levels  are  derived  which  are  valid  after  the  initial  time  transients  have 
disapeared.  Direct  numerical  simulations  of  two-dimensional  isotropic  compressible 
turbulence  are  performed  to  validate  the  results. 


II.  Theoretical  Analysis 


A.  Problem  Formulation 


The  compressible  Navier  Stokes  equations  are,  in  non-dimensional  variables, 


!  +  V.(,u)  =  0, 


d(pu) 

dt 


+  V-(puu)  = 


dP 

—  +  u  ■  VP  +  7PV  U  = 
ot 
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where 


a  -  2fj, 


Li(Vu  +  VuT)-i(Vu)/ 


(2) 


is  the  viscous  stress  tensor  (with  bulk  viscosity  assumed  zero),  I  is  the  identity  tensor, 
and 

$  =  I(Vu  +  VuT)  :  a  (3) 

Z 

is  the  viscous  dissipation  function.  The  equation  of  state  is 

P  =  pT.  (4) 


The  density,  velocity,  temperature  and  pressure  are  respectively 

PR 

Ur’ 

T* 

Tr’ 

P* 

RgpR.TR 

where  the  dimensional  reference  values  are  denoted  by  the  subscript  R  and  Rg  is  the 
universal  gas  constant.  Dimensional  variables  have  an  asterisk  superscript.  Distance 
and  time  are  scaled  respectively  with  respect  to  Lr  and  tR  —  Lr/Ur.  In  the  text,  x 
refers  to  the  cartesian  position  vector. 

The  Prandtl  number  Pr  —  (fiRCp)/ kr,  the  Reynolds  number  Re  =  (prUrLr)/ hr, 
7  is  the  ratio  of  specific  heats.  Viscosity  and  conductivity  are  denoted  by  /x  and  k 
respectively  and  are  scaled  with  respect  to  hr  and  kr. 


P  = 

u  = 

T  = 
P  = 


The  reference  Mach  number 


is  related  to  the  time-dependent  turbulent  Mach  number 


(5) 


Mt 


7  RgT* 
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Mr\  < 


(6) 


The  objectives  of  the  following  analysis  is  first  to  classify  the  various  types  of 
equilibrium  turbulent  regimes  (distinguished  by  the  presence  or  absence  of  shocks  and 


also  by  the  fraction  of  the  total  kinetic  energy  solely  due  to  compressibility  effects), 
and  second,  to  predict  the  range  of  initial  conditions  that  leads  to  each  one.  The 
effect  of  viscosity  is  felt  either  on  a  viscous  time  scale  (much  greater  than  the  acoustic 
time  scale),  or  during  the  formation  of  shocks.  In  the  latter  case,  viscosity  serves  to 
prevent  the  formation  of  singularities.  Although  there  is  a  distinct  possibility  that 
shocks  will  form  as  a  result  of  certain  types  of  initial  conditions,  viscosity  does  not  help 
initiate  the  processes  (wave  steepening)  which  eventually  lead  to  shocks.  Rather,  the 
viscosity  diffuses  the  sharp  gradients  to  generate  shocks  of  finite  thickness.  The  above 
considerations  lead  us  to  neglect  all  viscous  effects  in  the  theoretical  formulation  of 
the  full  viscous,  Navier-Stokes  equations.  This  assumption  will  be  verified  a  posteriori 
by  the  direct  numerical  simulations. 


In  the  following  analysis,  the  turbulent  Mach  number  is  assumed  to  be  very  much 
less  than  unity,  so  that  the  sound  velocity  is  much  greater  than  the  flow  velocity. 
Under  these  circumstances,  the  acoustic  time  scale  is  much  less  than  the  convective 
time  scale,  which  in  turn  is  much  less  than  the  viscous  time  scale  (if  the  Reynolds 
number  is  sufficiently  large).  As  will  be  shown  later,  the  different  regimes  emerge  on 
the  acoustic  time  scale,  during  which  time  any  inconsistency  in  the  acoustic  compo¬ 
nent  of  the  flow  is  washed  away.  In  other  words,  time  boundary  layers  only  extend 
over  a  time  period  of  O(Mr). 


Dropping  the  viscous  and  heat  conduction  terms,  Eqs.  (1)  become  (in  nonconser¬ 
vative  form) 

dp 

-57  +  u  •  Vp  +  pV  u  =  0, 
at 

^+uVu)  +  ^vp  =  0'  (7) 
dP 

— — b  u  •  VP  =  — qPV-u. 
at 

The  initial  conditions  are  arbitrary  with  the  restriction  that  both  Mr  and  the  root 
mean  square  (rms)  pressure,  ||p(x,  0)||,  be  much  less  than  unity.  For  the  remainder  of 
this  analysis,  we  assume  that  for  homogeneous  turbulent  flows,  the  maximum  norm 
and  the  L2  norm  are  of  the  same  order  of  magnitude. 


B.  Asymptotic  Analysis 

The  system  of  equations  (7)  is  now  investigated  by  an  asymptotic  analysis  which 
assumes  Mr  <<  1.  The  pressure  P  is  decomposed  into  a  mean  and  a  fluctuating 
component: 

P  =  1  +  P,  (8) 

where 

I  Ml  «  1.  (9) 
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For  simplicity,  we  neglect  any  consideration  of  the  entropy  mode  and  assume  the 
homentropic  condition 

P  =  p\  (10) 

The  density  is  expanded  in  powers  of  p  according  to 

P  =  1  +  7_1P+  0(|!p||2).  (11) 


Using  the  expansions  (8)  for  pressure  and  (11)  for  density,  Eqs.  (7)  become  (to  first 
order  in  ||p||) 


<9 u  1 

W  +  u'Vu  +  WRVp  =  °' 


dp 

+  u  •  Vp  -fi  7V-u  =  0. 
at 


(12) 


At  t  —  0  the  initial  data  is 


u(x,0)  =  ii0(x), 
p(x,0)  =  po(x). 


(13) 


It  is  convenient  for  the  purpose  of  analysis  to  split  the  velocity  vector  according 
to 

u  =  u7  +  uc,  (14) 

where  the  solenoidal  component  u7(x,t)  satisfies  the  time-dependent  incompressible 
Navier-Stokes  equations 


3u7 

~dt 


+  u7  •  Vu7 


Vp7, 


Vu7  =  0, 

with  initial  conditions 

u7(x,0)  =  u7(x), 

Vu7  =  0. 

The  incompressible  pressure  p7  satisfies  the  Poisson  equation 

V2p7  =  -Vu7  :  Vu7 

from  which  it  is  inferred  that 

p7  =  0(u72). 


(15) 


(16) 


(17) 

(18) 
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The  initial  solenoidal  velocity  is  determined  from  the  decomposition 


u0(x)  =  Uo(x)  +  u£(x).  (19) 

If  the  flow  is  homogeneous,  it  is  easy  to  show  that  Uq(x)  and  Uq(x)  are  unique 
once  u0(x)  is  specified.  When  the  flow  is  non-homogeneous,  an  arbitrary  potential 
function  can  be  added  to  the  decomposition  (19).  The  transverse  velocity  component, 
u7  is  divergence- free  for  all  time  (since  it  satisfies  the  time-dependent  incompressible 
Navier-Stokes  equations).  Although  uc  is  initially  vorticity-free,  it  acquires  a  small 
degree  of  vorticity  subsequently.  Kreiss,  Lorenz  and  Naughton  (1990)  have  estimated 
this  vorticity  to  be  O(M^). 

Substitution  of  Eq.  (14)  into  Eqs.  (12)  and  manipulation  of  the  resulting  expres¬ 
sions  using  Eqs.  (15)  yield 

^ — h  u7  •  Vuc  +  uc  ■  Vu7  +  uc  •  Vuc 
ot 

^  +  u7  •  Vp  +  uc  •  Vp  +  7V-uc 
ot 

The  momentum  equation  in  Eqs.  (20)  is  simplified  if  the  perturbation  pressure  is 
decomposed  according  to 

p  =  7  Af£p7  +  8pc.  (21) 

This  particular  decomposition  removes  p1  from  the  evolution  equation  for  uc .  At 
t  =  0, 

pJ(x,0)  =  pI0(x)  (22) 

is  uniquely  determined  from  the  incompressible  evolution  equations,  so  that  there 
is  no  ambiguity  in  the  computation  of  pc(x,0).  The  uniqueness  of  the  pressure 
decomposition  is  maintained  for  all  time.  The  parameter  8  is  defined  to  insure  that 

||pc(x,0)|j  =  l.  (23) 

The  total  initial  velocity  u  is  also  normalized  to  unity,  i.e.  ||u0(x)||  =  1.  This  is  most 
easily  accomplished  by  choosing 


(VP-7MJV/), 


(20) 


=  0. 


Ur  =  IK(x)||.  (24) 

Therefore,  both  Uq  and  Uq  are  0(1).  From  Eq.  (18),  | |pol |  =  0(1).  Incompressible 
variables  vary  on  a  slower  time  scale  than  the  acoustic  time  scale,  which  is  the  time 
scale  of  interest,  so  that  their  order  of  magnitude  will  remain  invariant  in  time.  In 
contrast,  we  will  demonstrate  that  certain  combinations  of  initial  conditions  can  lead 
to  time  boundary  layers  which  can  change  the  order  of  magnitudes  of  compressible 
variables  on  a  time  scale  of  O(Mr). 
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Substitution  of  Eqs.  (14),  (21)  into  the  system  (12)  leads  to  the  evolution  equations 

^  +  u1  ■  Vuc  +  uc  •  Vuc  +  Auc  +  -4tt^PC  =  0. 
dt  7  Mft 

+  u1  •  Vpc  +  uc  •  Vpc  +  Buc  +  1 V  uc  =  G, 
at  o 

for  the  compressible  components  of  the  flow  variables,  where 

Amc  —  uc  •  VuJ, 

Buc  =  •  Vp', 

G  =  + 

For  reference,  the  orders  of  magnitude  for  A ,  B  and  G  at  t  =  0  are 

A  =  0(1), 

B  =  O(^), 

G  =  O0?&). 

0 

Since  A,  B,  G  only  depend  on  incompressible  quantities,  their  order  of  magnitude  is 
constant  over  the  time  scale  of  interest.  The  initial  data  is 

uc(x,0)  =  u0(x)  Uq(x)  —  u^(x), 

PC(x,  0)  =  Po(x)  -  =  Po(x)- 

Note  that  Uq  is  vorticity-free.  The  pressure  P(x)  is  now  given  by 

P(x,  t)  =  1  +  7 M2Rpl{x,  t )  +  6pc(x ,  t).  (25) 

It  must  be  noted  that  although  pc  and  uc  are  0(1)  at  t  =  0,  there  is  no  guaranty 
that  they  will  remain  so.  Should  they  exceed  that  order,  they  will  be  scaled  down, 
otherwise  terms  might  be  neglected  in  the  governing  equations  which  should  be  kept. 


C.  Regime  Classification 


Freezing  coefficients  in  Eqs.  (25)  makes  it  clear  that  there  are  two  different  time  scales 
present,  i.e.  the  convection  scale  of  order  0(1),  and  the  sound  wave  scale  of  order 
O(Mr).  Since  the  initial  compressible  velocity  field  u£ (x)  is  vorticity-free,  we  expect 
uc  to  vary  on  the  fast  acoustic  time  scale.  Therefore  we  introduce  a  new  time  variable 
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0, 


(27) 


In  terms  of  t',  Eqs.  (25)  become 

^  +  MR[(uI  +  uc)-Vuc  +  Auc]  +  jL-Vpc  = 

+  Mr  [(u7  +  uc)  •  Vpc  +  Buc]  +  -jp-V-uc  =  MrG. 

The  system  of  equations  (27)  is  hyperbolic  and  for  8  <<  Mr  or  8  »  Mr  it  is  highly 
asymmetric.  Therefore  the  solution  will  in  general  grow  rapidly,  thus  invalidating 
the  original  scalings  (see  appendix).  Only  with  the  symmetrized  version  of  Eqs.  (27) 
can  we  be  sure  that  the  solution  at  t  >  0  is  bounded  by  the  order  of  magnitude  of 
the  solution  at  t  =  0.  Therefore  we  symmetrize  the  system  in  such  a  way  that  the 
initial  data  are  not  “scaled  up”  because  otherwise  nonlinear  terms  which  are  initially 
negligible  may  not  remain  so.  The  change  of  variables  which  leads  to  symmetry  is 
dependent  on  whether  the  coefficient  of  Vpc  in  Eq.  (27)  is  less  or  greater  than  unity. 
In  terms  of  8,  we  therefore  distinguish  two  cases. 

8  <  Mr') 


This  regime  is  characterized  by  low  levels  of  initial  acoustic  pressure  fluctuation. 
We  first  introduce  a  new  dependent  variable 


(28) 


which  insures  that 
Eqs.  (27)  become 


I  Ip? 'II  <  I  boll  =  0(1)- 


d\xc 

~W 


+  Mr[( u1  4-  uc)  •  Vuc  +  Auc]  +  Vpc' 


dt 

with  initial  data 


~  +  Mr\(  u'  +  uc)  ■  Vpc  +  zjnBuC  1  +  v'“c 


o, 


(29) 


(30) 


uc(x,0) 

Pc'(x,0) 


«o(x), 


By  construction,  uc(x,  0)  and  pc'(x,0)  are  0(1).  If  they  should  be  much  less  than 
unity,  the  dependent  variables  are  rescaled  so  that  they  become  0(1),  in  which 
case  the  coefficients  of  Eqs.  (30)  are  replaced  by  smaller  coefficients,  and  the  ap¬ 
proximations  below  become  more  valid.  (The  notation  /  =  O(M)  is  equivalent  to 
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liniAf_o  f{M)/M  ^  0  <  oo.  On  the  other  hand,  the  0{M)  symbol  only  states  that 
the  limit  is  finite,  possibly  zero.)  The  system  (30)  is  a  well-behaved  system  of  dif¬ 
ferential  equations.  However,  one  must  check  whether  the  convective  terms  remain 
negligible  when  t'  =  0(1).  It  is  easily  shown  that  on  the  fast  time  scale,  the  lowest 
order  evolution  equations  are  given  by  the  linear  wave  equations 


duc 

~w 


+  vPc 


0, 


with  the  initial  data 


dpc> 

~W 


+  Vuc 


0, 


uc(x,0)  =  u£(x), 
pc'(x’0)  = 

On  the  convective  time  scale,  Kreiss  et  al.  (1990)  have  shown  that  the  wave  equation 
still  describes  the  acoustic  component  of  the  flow. 


From  the  wave  equation,  it  is  now  easy  to  predict  the  order  of  magnitude  of  uc 
and  pc  as  a  function  of  initial  magnitudes.  After  a  time  t'  =  0(1),  both  uc  and  pc' 
will  have  the  same  order  of  magnitude,  given  by 

uc(x,*')  »  pc'(x,t')  =  0(max(u£(x),p£'(x))), 

=  0(max(ug(x), 

The  results  are  summarized  in  table  1  for  all  different  initial  estimates  of  velocity. 


PC  0 

PCo o 

O(l) 

0(M’R) 

0(1) 

0(1) 

0(1) 

0(max(M£,  ,£j) 

o(®n 

0(max(^-f-,  1)) 

Table  1:  Order  of  magnitude  of  ~quilibrium  levels  of  compressible  pressure  and 
acoustic  velocity  as  a  function  of  initial  levels. 


With  the  chosen  initial  normalization  for  the  velocity,  u;  must  be  0(1),  unless 
uc  =  0(1).  The  choice  of  the  initial  magnitude  of  the  incompressible  velocity  compo¬ 
nent  partly  determines  the  initial  variation  of  the  ratio  of  compressible  kinetic  energy 
to  total  kinetic  energy. 

8  >  Mr  7 
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Now  we  introduce  the  new  variable 


C>  Mr^  g 

u°  =  — —  u° 


which  insures  that 

IK'II  <  lluoll  =  0(1). 

Using  Eq.  (31),  Eqs.  (27)  become 


+  {Mru!  +  -uc')  •  Vuc'  +  MrAvlc'  +  Vpc  =  0, 


??-  +  (Mru1  +  -uc')  •  Vpc  +  -Buc'  +  V-uc'  =  MrG 

ot'  7  7 


with  the  initial  data 

uc'(x,0)  =  u?(x), 

(34) 

PC(x,  0)  =  p£(x)  =  0(l). 

Once  ’’gain,  if  the  initial  data  in  Eq.  (34)  should  be  much  less  than  unity,  a  rescaling 
of  velocity  and  pressure  would  simply  decrease  the  coefficients  of  the  quadratic  terms. 
As  long  as  uc,(x,f')  =  0(1),  Eqs.  (33)  reduce  to 

^-+-uc'-Vu  c'  +  Vpc  =  0, 

at'  7 

(35) 

^  +  -uc'-Vpc  +  -Buc'  +  V-uc'  =  0, 

at1  7  7 

with  initial  conditions  given  by  Eq.  (34).  When  t'  =  0(1), 

lluC'(x,0ll  =  ilPC(*,f')ll. 

=  0(max(p£,u£')),  (36) 

=  0(1). 

Equations  (31)  and  (36)  imply  that 

uc(x,t')  =  (37) 


When  S  —  0(1),  the  convective  terms  balance  the  time  derivative  terms,  in  which 
case  wave  steepening  may  occur  on  a  t1  =  0(1)  time  scale,  and  there  is  a  propensity 
for  shocks  to  form.  Note  that  when  shocks  do  occur,  the  turbulent  Mach  number,  Mt 
is  no  longer  small,  but  has  increased  by  a  factor  1  /Mr  from  its  initial  value. 
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D.  Results 


In  the  previous  section,  we  established  the  conditions  for  two  separate  asymptotic 
regimes,  depending  on  whether  the  ratio  of  8  to  7 Mr  was  lesser  or  greater  than  unity. 
We  now  wish  to  deduce  some  physical  consequences  from  the  results  of  the  previous 
section.  To  this  end,  we  consider  various  levels  of  initial  rms  values  for  uc0  and  8. 
The  rms  of  u  Ms  fixed  at  unity.  Only  when  j|ug  ||  =  0(1)  can  ||u7||  <  0(1).  It  is  then 
easily  shown  that  the  rms  level  of  u7  can  be  increased  to  unity  without  affecting  the 
order  of  magnitude  estimates  at  t'  =  1. 


Table  2  summarizes  the  theoretical  findings  from  the  previous  section.  To  properly 
interpret  the  table,  replace  a  cell  value  of,  for  example,  3  by  0(Mr).  A  subscript  Q,*, 
refers  to  the  value  of  a  variable  after  the  initial  transients  have  died  away  {t'  =  0(1)), 
also  referred  to  as  an  equilibrium  state.  When  viscous  effects  are  taken  into  account, 
the  equilibrium  state  of  all  quantities  varies  over  the  viscous  time  scale.  They  also 
fluctuate  in  time  due  to  higher  order  convective  effects  that  have  been  neglected  to 
lowest  order  in  the  current  expansion.  For  example,  on  the  incompressible  convective 
time  scale,  the  terms  Mr u7  •  Vuc  induce  small  fluctuations  of  the  variables  about 
the  equilibrium  state  of  the  system. 


The  rms  levels  of  Uq  and  6  are  given  in  the  first  two  columns  of  table  2.  The  fifth 
and  sixth  columns  contain  the  new  quantity 


X  = 


u 


C112 


u  I!2  +  !lu 


I\\2 


(38) 


which  is  the  fraction  of  the  total  kinetic  energy  solely  due  to  acoustic  effects.  Note 
that  at  t  =  0,  the  denominator  is  unity  by  construction,  so  that 


lluoll  =  \A  -  Xo, 

and 

IWII  = 


(39) 

(40) 


From  Eq.  (37), 


Xoo 


1 

1  +  M£||u7||2’ 


(41) 


and  since  ||u7||  =  0(1),  8  =  0(1)  always  implies  that  Xoo  =  0(1).  Consequently, 
the  lower  the  initial  level  of  compressibility  at  a  fixed  Mr  (i.e.  the  lower  Xo),  the 
stronger  the  initial  transient  (which  extends  over  the  time  period  O(Mr)).  As  this 
imbalance  intensifies,  the  propensity  for  shocks  to  form  also  increases.  Table  2  re¬ 
flects  these  statements  in  the  rows  corresponding  to  8  =  0(1).  When  u£  =  0(1), 
u7  =  O(Mr)  ,  n  >  1.  However,  the  value  of  n  does  not  change  the  previous  conclu¬ 
sions.  Based  on  the  results  in  Sarkar,  Erlebacher,  Hussaini  and  Kreiss  (1989),  it  is 
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straightforward  to  show  that 


where 


x</l 

Aoo 


1  1  +(2 

2  1-0.5xo(1-£2)’ 


__  HpcH 

Mt7  ' 


(42) 


(43) 


u  o 

6 

- 

U  OO 

PC  OO 

Xo 

Xoo 

F0 

0 

0 

-1 

0 

0 

0 

2 

0 

1 

0 

0 

0 

0 

0 

0 

2 

0 

-1 

0 

0 

-2 

0 

3 

0 

-2 

0 

0 

-4 

1 

0 

-1 

0 

2 

0 

4 

1 

1 

0 

0 

2 

0 

2 

1 

2 

1 

0 

2 

2 

0 

1 

3 

1 

-1 

2 

2 

-2 

1 

4 

1 

-2 

2 

2 

-4 

2 

0 

-1 

0 

4 

0 

6 

2 

1 

0 

0 

4 

0 

4 

2 

2 

1 

0 

4 

2 

2 

2 

3 

2 

0 

4 

4 

0 

2 

4 

2 

-1 

4 

4 

-2 

2 

5 

2 

-2 

4 

4 

-4 

Table  2:  Estimates  of  equilibrium  levels  of  compressible  pressure  ( pc ),  (uc), 

X  is  a  function  of  a  wide  range  of  initial  levels  for  pc  and  uc.  Table  entries 
should  be  interpreted  as  0(MlI^ble  entrv). 

Two  separate  regimes  (different  from  the  mathematical  regimes  of  the  previous 
section)  are  apparent  from  table  2.  For  large  values  of  6,  ||p£j|  =  0(1),  while  uc 
increases  sharply.  This  sudden  increase  over  the  acoustic  time  scale  is  most  important 
when  6  is  large  and  Uq  is  small.  The  fact  that  this  increase  in  the  acoustic  velocity  is 
not  contingent  on  8  —  0(1)  indicates  that  an  initial  imbalance  in  compressible  energy 
can  occur  without  the  need  to  appeal  to  wave  steepening.  When  6  is  sufficiently 
small,  uc  remains  at  its  initial  level  and  large  pressure  waves  appear  in  the  flow. 
The  compressible  pressure  imbalance  becomes  more  severe  with  decreasing  and 
decreasing  6.  Equilibrium  initial  conditions  are  characterized  by  the  juncture  between 
the  two  aforementioned  regimes.  This  corresponds  to 

<  =  «4) 
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Using  Eq.  (39),  Eq.  (44)  becomes 


^  =  6(1).  (45) 

The  function 

m  =  (46) 

was  introduced  by  Sarkar  et  al.  (1989)  and  shown  to  play  a  fundamental  role  in  the 
description  of  compressible  turbulence.  The  last  column  in  table  2  tabulates  the  value 
of  Fo  =  F( 0)  for  the  various  cases.  It  is  now  easy  to  see  that  a  necessary  condition 
for  the  appearance  of  shocks  is  6  —  1  (see  previous  section).  If  the  latter  condition 
is  not  satisfied,  the  level  of  compressibility  will  remain  small  for  all  time,  which  is 
inconsistent  with  the  appearance  of  an  isotropic  distribution  of  shocks. 

It  is  interesting  to  establish  the  conditions  under  which  pc  is  a  good  approximation 
to  the  total  perturbation  pressure  p.  To  this  end  consider  the  relation  between  pc , 
p1  and  p: 

^iipciij  =  iipip  +  72mjiip'ii2  (47) 

which  is  justified  if  p(x)  are  p;(x)  are  assumed  to  be  decorrelated.  This  is  true  by 
construction  at  t  —  0  in  our  direct  numerical  simulation  code  (discussed  in  a  later 
section).  From  Eq.  (18),  ||p;||  and  ||u7||  are  related  by 

lip'll  =  a1||u/||2, 

=  «i(l  ~  Xo),  (48) 

where  ax  is  an  assumed  0(  1)  constant.  Using  Eq.  (48),  Eq.  (47)  becomes 

£2||PCH2  =  INI2  +  72M£a2(l  -  Xo)2-  (49) 

The  incompressible  component  of  the  pressure  field  is  therefore  negligible  when 

||po|]  >>  jM^a^l  -  xo)-  (50) 

This  result  will  be  verified  with  the  help  of  direct  numerical  simulations.  Equa¬ 
tion  (50)  tells  us  that  one  cannot  neglect  the  influence  of  pl  at  higher  levels  of  Mr 
and  at  lower  levels  of  the  initial  compressibility  ratio,  Xo- 

III.  Numerical  Simulation 

A.  Numerical  Algorithm 

Our  direct  simulations  of  two-dimensional  compressible  turbulence  are  based  on  Eqs.  ( 1 ) 
The  flow  is  assumed  to  be  isotropic.  Hence  a  double  Fourier  representation  is  appropri¬ 
ate.  Spectral  methods,  by  nature  of  their  high  accuracy  and  low  dispersive  and  dissi¬ 
pative  errors  are  ideally  suited  for  this  problem.  Spatial  derivatives  are  approximated 


by  Fourier  collocation  (Canuto,  Hussaini,  Quarteroni  and  Zang  1988).  In  each  coor¬ 
dinate  direction,  N  grid  points  are  used;  for  example  Xj  =  2-Kj/ N,  j  =  0,1,  -  ■  ■ ,  N  —  1 
in  the  x  direction.  The  derivative  of  a  function  .F(x)  with  respect  to  x  is  calculated 
from  the  analytic  derivative  of  the  trigonometric  interpolant  of  .P(x)  in  the  direction 
x.  Many  simulations  of  incompressible,  homogeneous  turbulence  have  used  a  Fourier- 
Galerkin  method.  The  compressible  equations,  however,  contain  cubic  rather  than 
quadratic  nonlinearities  and  true  Galerkin  methods  are  more  expensive  (compared 
with  collocation  methods)  than  they  are  for  incompressible  flow.  The  essential  dif¬ 
ference  between  collocation  and  Galerkin  methods  is  that  the  former  are  subject  to 
both  truncation  and  aliasing  errors,  whereas  the  latter  have  only  truncation  errors. 
As  discussed  extensively  by  Canuto,  et  al.  (1988),  the  aliasing  terms  are  not  signifi¬ 
cant  for  a  well-resolved  flow.  However,  care  is  needed  to  pose  the  relevant  equations 
for  a  collocation  method  in  a  form  which  ensures  numerical  stability.  For  this  reason, 
the  convective  term  in  the  momentum  equations  is  actually  used  in  the  equivalent 
form 

-[V-(/>uu)  +  pu  •  Vu  +  uV-(pu)].  (51) 

As  noted  by  Feiereisen,  Reynolds  and  Ferziger  (1981),  this  form,  together  with  a 
symmetric  differencing  method  in  space  (for  example  Fourier  collocation),  ensure 
conservation  of  mass  and  momentum.  Moreover,  energy  is  conserved  for  the  semi¬ 
discrete  inviscid  equations. 


We  are  interested  in  the  low  subsonic  regime  Mr  <<  1.  The  sound  speed  is  in 
this  case  much  greater  than  the  flow  velocity  and  this  imposes  a  severe  restriction 
on  the  time  step  of  any  explicit  time  marching  numerical  scheme.  To  remove  this 
constraint,  a  splitting  method  is  adopted.  The  first  step  integrates  implicitly 


dp 


=  0, 


dt 

-Jp  +  ^[V-(^uu)  4-  u  •  V(pu)  +  pu  •  Vu]  =  ^Ver,  (52) 

dP 

—  +  u  •  VP  -)-  7PV-u  -  CqV’(^u)  = 


7 


RePr 


v2r 


■  7(7~1  )M2r 
Re 


from  which  the  acoustic  effects  that  vary  on  the  fast  time  scale  0{MR)  have  been 
removed.  The  second  step  integrates  the  “acoustic  equations”  which  vary  on  the  fast 
time  scale.  These  equations  are 


+ v(pu) 


dpu 

~di 

dP 

dt 


+ 


1 


VP 


+  coV  ■  {pu) 


=  0, 

=  0, 

-  0. 


(53) 
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where  cq  is  the  root  mean  square  (rms)  of  the  sound  speed. 


This  splitting  is  employed  at  each  stage  of  a  low-storage  third-order  Runge-Kutta 
method.  To  minimize  the  errors  due  to  the  large  implicit  time  step,  Eqs.  (53)  are 
integrated  analytically.  In  Fourier  space  the  system  (53)  become 


dp 

5F  +  lkm  =  °- 
9m  tk  « 

~dt+^M%  =  °’ 

dP  2.  -  n 

—  +  ic0k  •  m  =  0, 


(54) 


where  m  =  pu  and  Fourier  transformed  quantities,  which  depend  upon  the  wavenum¬ 
ber  k,  are  denoted  by  a  circumflex.  Equations  (54)  are  solved  exactly  in  Erlebacher, 
Hussaini,  Speziale  and  Zang  (1987)  and  are  rewritten  here  for  completeness: 


p(’)  = 


i  ^  r  A  /CokAt,.  *  .  .  C0kAt,  .  ■j. 

’ +  4^  cos  W’ +  m  W’  L 


-  (2) 
rrr  ' 


i<2) 


-  (1) 

—  m'  '  — 


zk 


c0kMRs/i 


'  = 


2  .Co kAt,.  *  .  .CokAt,. 
A  co<in—i=)  +  B  sm(77 — ;=)> 


r  j  «  /  Cok&ta  .  * 

iA  sm(77-^)  -  B  cos(tt—7=)  +  Bl 


'MRyfj* 


(55) 


MR^yJ 


Mr^ y' 


where  k  =  |k|  is  the  magnitude  of  the  Fourier  wavevector,  A  =  P^\  B  =  i^-k  •  m^. 
Superscripts  1  and  2  relate  to  the  end  of  the  first  and  second  fractional  step  respec¬ 
tively.  The  effective  time-step  of  the  Runge-Kutta  stage  is  denoted  by  At,.  The 
advantage  of  this  splitting  is  that  the  principal  terms  responsible  for  the  acoustic 
waves  have  been  isolated.  Since  they  are  treated  semi-implicitly,  one  expects  the 
explicit  time-step  limitation  to  depend  upon  v  +  |c  -  c0|  rather  than  v  +  c.  (Although 
there  is  also  a  viscous  stability  limit  for  the  first  fractional  step,  it  is  well  below  the 
advection  limit  in  the  cases  of  interest.)  This  is  clearly  efficient  at  low  Mach  num¬ 
bers.  Since  the  second  fractional  step  is  integrated  analytically,  it  does  not  generate 
dispersion  or  dissipation  errors.  The  only  errors  incurred  are  time  splitting  errors. 


If  one  is  truly  interested  in  the  detailed  time-evolution  of  the  acoustic  components, 
the  time-step  must  be  small  enough  to  resolve  the  temporal  evolution  of  these  waves. 
Note  however,  that  because  the  acoustics  are  treated  exactly  in  the  second  step,  the 
large  time  step  simulation  still  produces  accurate  results,  but  the  data  is  not  available 
at  intermediate  times. 


The  expected  stability  limit  of  this  two-dimensional  Fourier  collocation  method 
for  the  compressible  Navier-Stokes  equations  has  the  form 


At  <  a  min( 


Ax 


+ 


Ay 


u|  +  |c-c0|  |u| -f  |c  -  c0| 


)• 


(56) 
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For  the  third-order  Runge-Kutta  method  employed  here  a  —  0.54.  For  accuracy 
purposes,  the  results  herein  are  all  based  on  a  =  0.1. 


B.  Initial  Conditions 

The  initial  conditions  for  the  direct  numerical  simulations  presented  here  are  similar 
to  those  presented  by  Passot  and  Pouquet  (1987).  They  are  sufficiently  general  to 
produce  the  various  turbulent  regimes  we  expect  to  occur  when  Mr  <<  1.  One  mus* 
separately  specify  the  reference  Mach  number  Mr,  the  Reynolds  number,  Re,  the 
Prandtl  number,  Pr,  the  rms  levels  of  u£ ,  u!0,  j|p0||,  ||T0||,  and  the  autocorrelation 
spectrum  for  p,  T,  uc  and  u7.  We  now  detail  the  procedure  used  in  the  numerical 
code. 

We  first  generate  random  fields  in  physical  space  in  the  range  [—0.5, 0.5]  for  u0) 
po,  and  To.  Next,  the  velocity  is  decomposed  into  the  sum 

u0(x)  =  u£(x)  +  u7(x)  (57) 

where 


VuJ  =  0, 

Vxuq  -  0. 

This  decomposition  is  unique  for  homogeneous  flows  and  is  accomplished  in  Fourier 
space  according  to  the  prescription 


After  transforming  the  velocity  components  back  to  physical  space,  we  rescale  uc, 
u7,  p,  and  T  to  impose  a  prescribed  autocorrelation  spectrum.  The  autocorrelation 
spectrum  £pp(k)  of  a  variable  p(x)  is  defined  as 

J  p(x)p(x)dx  =  J  Epp(k)dk.  (58) 

For  example,  the  2D  energy  spectrum  for  u,  i.e.  E(k),  is  the  autocorrelation  spectrum 
l?u-u- 

We  now  impose  the  spectrum 

_  k* 

E(k)  =  k4e  (59) 

on  each  of  the  variables.  The  procedure  for  imposing  the  required  spectrum  E(k)  is 
explained  using  the  density  as  an  example.  One  scans  the  wavenumber  space,  and 
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for  each  pair  ( kx ,  kv),  the  corresponding  spherical  shell  is  determined.  The  ith  shell  is 
defined  by 

(i-^<k<i+1-).  (60) 

where  k  =  <Jk2x-\-  k2.  All  spherical  shells  have  unit  width.  The  0i/l  shell  is  empty 
because  k  =  (0,0)  is  the  mean  component  which  does  not  contribute  to  the  autocor¬ 
relation  spectrum.  The  contribution  of  the  ith  shell  to  the  autocorrelation  spectrum 
thus  becomes 

K  =  E  !«k)l2  (61) 

*-  j<*<»+ j 

where  E*  is  the  discrete  representation  of  E*(k)dk.  (For  2-D  flows,  spherical 

should  be  interpreted  as  circular.)  Defining  E{  as 


Ex  =  E(z) 

where  E{i)  is  obtained  from  Eq.  (59),  the  density  is  rescaled  according  to 


K')  < —  K') 


(62) 


(63) 


in  order  to  impose  the  desired  spectrum.  In  a  similar  way,  the  autocorrelation  spectra 
for  uc,  u 7  and  T  are  calculated.  For  the  velocity,  the  procedure  is  identical  to  that 
described  for  p,  except  that  the  scalar  p  is  replaced  by  a  vector,  and  |p(k)|2  is  replaced 
by  |u(k)|2. 

It  is  still  possible  to  rescale  these  quantities  by  an  arbitrary  constant  in  physical 
space  without  changing  the  properties  of  its  autocorrelation  spectrum.  Therefore,  wc 
impose  given  rms  levels  for  p,  T  and  u  by  an  appropriate  rescaling.  There  still  remains 
a  degree  of  freedom  on  the  velocity,  since  uc  and  u7  can  be  weighted  independently 
without  affecting  the  rms  of  u.  Therefore,  we  specify  the  initial  level  of  compressibility 
in  the  flow  according  to 

/  u °2dx 


A  /(uc2  +  u'!Vx'  {64) 

The  denominator  of  Eq.  (64)  is  the  total  kinetic  energy  of  the  system  since  /  uc  u7  =  0 
by  construction. 

The  mean  values  of  p  and  T  are  then  added  to  the  thermodynamic  variables  to 
obtain 


p  < —  1  +  p, 

T  < —  1  +T. 


The  reference  pressure  in  the  code  is  chosen  to  be  Pr  =  prUr\  consequently  the 
nondimensional  pressure  is  then  calculated  from 


P  = 


PT 


(65) 
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where  p  and  T  are  respectively  the  non-dimensional  density  and  temperature.  The 
turbulent  Mach  number  Mt  is  related  to  the  reference  Mach  number  Mr  according 
to  Eq.  (6). 


IV.  Numerical  Validation  of  Theory 


In  this  section,  we  show  results  of  2-D  direct  numerical  simulations  on  grids  of  64  x 
64  to  validate  the  theory  described  in  the  preceding  sections.  This  validation  is 
accomplished  by  comparing  the  value  of  x  after  several  acoustic  periods  as  predicted 
by  the  computation,  against  that  given  by  the  theory.  The  parameter  F,  defined  by 
Eq.  (46)  is  also  computed  and  compared  with  numerical  results. 

Given  the  large  parameter  space,  in  every  run,  the  density  and  temperature  fluctu¬ 
ation  rms  levels  are  set  equal  to  each  other.  An  extensive  parameter  range  is  covered 
by  considering  all  combinations  of 

Xo  =  (0.1, 0.5, 0.9), 

Mr  =  (0.03,0.08,0.3), 

IHI  =  \\T\\  =  (0.01,0.05,0.1). 

Cases  are  referred  to  by  a  3  digit  combination.  For  example,  case  123  refers  to 
Xo  =  0.1,  Mr  =  0.08  and  ||p||  =  ||T||  =  0.10.  Simulations  are  run  at  Re  =  150, 
||u||  =  1,  k0  =  10.  This  corresponds  to  a  microscale  Reynolds  number  R\  =  20.  The 
various  diagnostics  are  computed  at  every  iteration. 

The  quantities  considered  for  comparison  between  the  numerical  computations 
and  the  theory  results  are 

(66) 

and  x-  Both  these  quantities  are  computed  by  averaging  them  over  several  consecutive 
iterations:  x  is  averaged  over  iterations  5  to  10  (xs-io)  and  10  to  15  (xio-is),  while 
F  is  averaged  over  iterations  10-15  (F10-15).  Their  equilibrium  values  Xoo  and  Foo 
are  computed  by  averaging  their  values  over  the  last  60  iterations  of  the  400  iteration 
numerical  simulation.  The  two  sets  of  averages  for  x  give  an  idea  of  the  convergence 
rate  of  x  towards  the  theoretical  prediction  xth •  If  X  is  averaged  over  too  few  time 
steps,  the  value  could  be  either  over  or  under  predicted  since  x(0  is  an  oscillatory 
integral  which  tends  to  unity.  On  the  other  hand,  if  x  is  averaged  over  too  many 
steps,  the  effects  of  viscosity  will  be  noticeable,  and  agreement  with  x«  will  not  be 
reached.  Averages  over  five  iterations  prove  adequate. 

Table  3  summarizes  the  simulation  results.  Note  that  Xio-is  is  in  closer  agreement 
with  Xto  u  than  is  X5-10  On  the  other  hand  the  average  x  decays  as  a  function  of 
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X5-10 

XlO-15 

Xoo 

wan 

0.11 

0.08 

KgS 

0.62 

0.44 

ns 

0.86 

0.76 

0.06 

0.06 

0.05 

0.25 

0.21 

0.14 

0.58 

0.52 

0.32 

0.08 

0.05 

0.06 

0.09 

0.07 

0.06 

0.11 

0.13 

0.09 

0.24 

U.  (  ( 

U.  1  1 

0.59 

|  0.94 

0.92 

0.84 

0.30 

0.35 

0.21 

0.47 

0.47 

0.29 

0.73 

0.69 

0.29 

0.40 

0.27 

0.47 

0.40 

0.29 

0.20 

0.42 

0.35 

0.23 

0.84 

0.84 

0.67 

0.95 

0.95 

0.88 

0.99 

0.98 

0.92 

0.78 

0.83 

0.63 

j  0.86 

0.86 

0.70 

0.94 

0.92 

0.81 

0.85 

0.76 

0.63 

0.85 

0.77 

0.63 

0.86 

0.80 

0.64 

*/NpoII 


1.00 

1.00 

1.00 

1.05 

1.00 

1.00 

5.10 
1.39 

1.10 


1.00 
1.00 
1.00 
1.01 
1.00 
1  00 
2.93 
1.12 
1.03 


0 
0 
0 
1.00 
1.00 
1.00 
1.12 
1.00 
1.00 
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time  in  part  due  to  viscous  effects  which  is  reflected  by  the  low  values  of  Xoo  compared 
to  Xx  in  fable. 


Sarkar  et  al.  (1989)  showed  analytically  that,  for  low  Mach  number  turbulence, 
the  parameter  F  approaches,  and  then  oscillates  about  a  value  of  —  1.  Table 
3,  regarding  F <*,  are  in  agreement  with  this  result.  The  value  of  F  depends  on  the 
compressible  component  of  pressure,  6pc ,  not  p.  Nonetheless,  it  is  p,  not  pc  which  is 
used  to  compute  F  for  table  3.  Furthermore,  the  time  interval,  required  to  obtain  a 
good  average  over  the  acoustic  oscillations  also  increases  with  Mr.  This  explains  why 
■Fio-15  is  not  close  to  unity  at  the  higher  values  of  Mr.  The  last  column  measures 
the  departure  of  the  initial  compressible  rms  pressure  from  the  initial  rms  pressure. 
The  cases  where  6/||p0||  >  1  correlate  well  with  a  non-equilibrium  value  of  F  after  15 
iterations.  On  the  viscous  time  scale,  F  is  nonetheless  close  to  unity  in  all  cases,  which 
reflects  that  the  turbulent  Mach  number  has  decreased  substantially  (see  Eq.  (50)), 
i.e.  p1  has  become  negligible  with  respect  to  pc . 

The  ratio  of  Xio-is  to  is  furthest  from  unity  when  Mr  =  0.3,  which  is  explained 
by  the  longer  acoustic  time  scale.  Since  the  averages  are  computed  over  a  fixed  number 
of  iterations,  the  sample  length  at  Mr  =  0.3  is  not  adequate,  and  the  prediction  for 
Xoo  is  not  very  good.  The  lower  value  of  Xoo  with  respect  to  x^  arises  because  x 
decays  on  the  viscous  time-scale.  On  the  other  hand,  F  remains  close  to  unity,  and 
is  basically  unaffected  by  the  viscosity  terms. 

Turbulent  simulations  with  initial  flowfields  for  which  F0  =  1  preclude  initial  time 
boundary-layers  which  might  otherwise  impose  stringent  restrictions  on  time-stepping 
and  code  accuracy.  We  now  derive  an  implicit  relation  for  x  as  a  function  of  Mr  and 
||p||.  To  this  end,  we  begin  by  collecting  the  required  formulas.  Henceforth,  all 
quantities  refer  to  initial  conditions  that  are  consistent  with  an  equilibrium  (F  =  1). 
Based  on  a  general  velocity  decomposition 

u  =  a(uc  +  /3ur)  (67) 


where  a  is  such  that  ||u|| 
as 


=  1,  the  compressible  to  total  kinetic  energy  ratio  is  written 

=  llucH2  =  W  rfyn 

X  ||uc||2+/?2||u/||2  7 2M2r  1  ] 


and  the  rms  of  the  compressible  pressure  field  becomes 


I!pcI!j  =  llpll!  +  yM>'||2, 

=  iipii2 + yMja^iiu'ir, 

where  we  have  assumed  that  p  and  pr  are  decorrelated.  By  construction,  they  are 
decorrelated  at  t  =  0  in  the  numerical  code  from  Eq.  (21)  More  realistically  how¬ 
ever,  p1  and  pc  are  probably  decorrelated  since  the  incompressible  and  compressible 
components  of  the  flowfield  evolve  quasi-independently. 
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From  the  previous  set  of  equations  one  obtains 

||p||2  =  72M£  [x  ~  cc\M2r(\  -  x)2 


(69) 


Contour  plots  of  x  are  shown  in  Figure  1  (aq  =  1).  It  is  seen  that  for  moderate  to 
high  values  of  x,  the  relationship  between  ||p||  and  Mr  is  linear,  which  simply  restates 
that  ||pc||  and  j jp) j  are  approximately  equal.  At  low  values  of  x,  !!pi!  decreases,  and 
eventually  reaches  zero.  In  other  words,  in  a  turbulent  flow  in  acoustic  equilibrium 
( F  =  1),  high  rms  pressure  at  high  turbulent  Mach  numbers  imply  that  x  's  close 
to  unity.  Since  a  necessary  condition  for  the  generation  of  shocks  is  8  —  0(1),  it 
follows  that  when  shocks  are  present  in  a  turbulent  compressible  homogeneous  flow 
in  acoustic  equilibrium,  t.he  compressible  part  of  the  flow  must  dominate. 


V.  Conclusions 

In  this  paper,  an  asymptotic  theory  that  is  based  on  an  initial- value  formulation  of  the 
equations  of  motion  for  compressible  homogeneous  turbulent  flow  has  been  presented. 
The  expansion  parameter  is  the  turbulent  Mach  number,  which  is  assumed  to  be 
small.  A  key  feature  of  the  theory  is  the  decomposition  of  the  turbulent  velocity  field 
into  two  physically  meaningful  components.  The  first  is  solenoidal  and  satisfies  the 
time-dependent  incompressible  Navier-Stokes  equations.  This  uniquely  determines 
the  second  component  which  is  initially  irrotational,  but  acquires  a  small  amount  of 
vorticity  at  later  times.  One  consequence  of  this  decomposition,  is  the  simultaneous 
separation  of  the  pressure  field  into  incompressible  and  compressible  components. 
The  incompressible  component  is  proportional  to  M^. 

To  lowest  order,  the  asymptotic  theory  describes  the  acoustic  part  of  the  flow  by  a 
wave  equation  which  is  coupled  to  the  solenoidal  flow  through  the  initial  conditions. 
At  later  times,  stronger  coupling  will  occur  through  the  neglected  convective  and 
viscous  terms  in  the  momentum  equation.  Order  of  magnitude  estimates  permit  an 
apriori  estimate  of  the  compressibility  levels  of  the  flow  as  a  function  of  the  initial 
conditions.  It  is  found  that  the  initial  flow  can  evolve  towards  two  basically  different 
states  (in  the  absence  of  shocks).  If  Fq  <  1,  the  pressure  fluctuations  remain  at  their 
initial  levels,  while  the  velocity  fluctuations  increase  in  such  a  way  that  x  becomes 
0(1).  On  the  other  hand,  if  Fq  >  1,  large  pressure  waves  develop  over  a  time  of 
0(Mr),  and  the  rms  level  of  the  acoustic  velocity  remains  constant.  This  transient 
behavior  is  well  described  by  the  solution  to  the  wave  equation,  and  has  been  checked 
against  extensive  two-dimensional  direct  numerical  simulations.  Note  however  that 
the  theory  is  also  valid  for  three-dimensional  turbulent  flows  and  spot  checks  have 
been  conducted. 

It  has  also  been  established  that  a  necessary,  but  not  sufficient,  condition  for  the 
generation  of  shocks  is  that  the  flow  must  be  in  an  initial  state  of  disequilibrium  (i.e. 
F  <  1)  which  results  from  an  initial  pressure  level  which  is  0(1). 
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Several  issues  need  further  clarification.  First,  it  is  still  not  clear  if  there  are  more 
stringent  initial  conditions  which  could  guarantee  the  emergence  of  isotropic  shock 
distributions  in  two  or  three-dimensional  compressible  turbulent  flows.  Moreover, 
it  is  not  clear  how  the  theory  should  be  modified  if  the  reference  Mach  number  is 
no  longer  small.  Second,  the  current  paper  only  treated  the  lowest  order  equations 
which  resulted  uom  the  asymptotic  expansions.  To  understand  the  evolution  of  the 
turbulent  flow  on  the  convective  and  viscous  time  scales,  it  is  necessary  to  keep  both 
viscous  terms,  and  higher  order  terms  in  the  Mach  number  expansion.  These  higher 
order  terms  will  allow  an  explicit  computation  of  the  effects  of  the  incompressible 
eddies  on  the  acoustic  waves  and  vice  versa.  These  two  issues,  are  now  being  studied 
and  will  be  the  focus  of  future  publications. 
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Appendix  A 


I.  Appendix 


We  want  to  explain  the  behavior  of  hyperbolic  systems  which  are  far  from  symmetric. 
Let  us  first  consider  the  system  of  ordinary  differential  equations 

u‘  = Au  =  ( o ) U|  u  =  {ul)'  <A1) 

with  initial  data 

u(0)  =  u0.  (A2) 

The  matrix  A  is  antisymmetric  (A7  =  —  A)  so  the  eigenvalues  fj.v  (j  =  1, 2)  are  purely 
imaginary  and  the  corresponding  eigenvectors  y:  are  orthonormal  with  respect  to  the 
scalar  product 


<  u,  v  >  =  u^vi  +  u2v2,  |u|2  =<  u,  u  >,  (A3) 

where  an  overbar  denotes  complex  conjugate.  The  general  solution  of  Eq.  (Al)  is 
given  by 

u(t)  =  a1e^ty1  +  a2e^ty2  (A4) 

where 

<7i=<yi,u0>,  o2  =<  y2,  u0  >  .  (A5) 

From  the  orthogonality  of  the  eigenvectors, 

iu(0!2  =  M2  +  M2  =  |u(0)|2,  (A6) 

and  the  amplitude  of  the  solution  does  not  change.  This  is  true  for  any  antisymmetric 
matrix  which  can  be  shown  by  the  energy  method 

d 

=  F  <u'u> 

=  <  Ut,  U  >  -f  <  U,Ut  > 

=  <  Au,  u  >  +  <  u,  Au  > 

=  -  <  u,Au  >  +  <  u,Au  > 

=  0. 
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By  the  previous  result, 

I  s(()!2  =  MOI’  +  KMl2 

=  |i,(0)|2  +  |u2(0)|! 

=  |ti!(0)|2  +  e|u2(0)|2. 

Assume  now  that  u^O)  =  0(1)  and  u2(0)  =  O'1).  The  general  form  (A4)  of  the 
solution  tells  us  that  in  general  both  {^(f)  -  0(1)  and  u2(t)  —  0(1).  Therefore  we 
obtain  for  the  original  variables 

i  (t)  =  ui(t)  =  0(1) 

u2(t)  =  e_1/'u2(t)  =  0(e~1/2) 

which  shows  that  the  amplitude 

KOI  =  O(€-1/2)|u(0)|  (aio) 

increases  dramatically  for  sufficiently  small  e  though  the  eigenvalues  of  A  are  still 
purely  imaginary.  We  can  also  express  the  result  in  the  following  way.  Assume  that 
Ui(0)  =  0  and  u2( 0)  =  1.  Then  | u(i ) |3  =  e  and 

K0I  =  1»  (A11) 

i.e.  for  this  initial  data,  the  amplitude  does  not  increase.  If  we  now  make  a  small 
perturbation  in  Ui(0),  say  ui(0)  =  77,  then  the  perturbation  will  be  amplified  such  that 
Ui(t)  «  0(1  +  t~ll2ri)  and  u2(i)  ~  0(1  +  e-1^2?/).  Thus  if  77  >>  e1^2  the  perturbation 
destroys  the  solution. 

There  are  no  difficulties  in  generalizing  the  results.  Consider  the  system 

ut  +  Au  =  0,  (A12) 
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where  A  is  an  n  x  n  matrix  and  u  is  a  vector  with  n  components.  Assume  that  there 
is  a  diagonal  transformation 


D  =  diag(di,-  ■  ■  ,dn),  0  <  dj  <  1,  (A13) 

such  that  DAD~X  is  antisymmetric.  Then  the  amplitude  of  the  solution  can  grow  by 
a  factor  maxjdj1,  i.e. 

|u(f)|  =  (^(maxd”1)  |u(0)|.  (A14) 

i 

In  particular,  if  we  initialize  the  data  such  that 

lu(t)|  ~  lu(0)|,  (A15) 

then  we  can  find  perturbations  tj  which  grow  to  the  order  max;  d~lr ]. 

After  Fourier  transform,  a  hyperbolic  system 

ut  =  Aux  (A16) 

becomes 

ue  =  iuAu.  (A17) 

The  eigenvalues  of  no  A  are  purely  imaginary  and  often  there  is  a  diagonal  transfor¬ 

mation  D  such  that 

lu)DAD~ 1  =  antisymmetric  matrix.  (A18) 

Then  we  can  apply  our  arguments.  As  an  example,  we  consider  the  wave  equation, 
written  as  a  first  order  system 

Ut  =  (l  *0  )  Ul  (A19) 

The  Fourier  transform  is  of  the  form  (A7). 
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Acoustic  Equilibrium 
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Figure  1:  |lp(Mn,  x)||  as  a  function  of  MR  for  different  values  of  x-  The  contour  levels 
of  x  aT  correspond  to  the  condition  that  F  =  1. 
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